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IS  abstract 


The  purpose  of  this  report  is  to  develop  a  general  algorithm  for  solving  the 
class  of  nonlinear  programming  problems  that  have  linear  constraints.  The  constraints 
can  be  either  equations  or  inequalities  and  the  variables  can  be  free  or  non-negative. 
The  objective  function  is  assumed  to  be  continuously  differentiable.  The  algorithm 
is  an  "effective**  second-order  method  in  that  slow  convergence  is  eliminated  without 
requiring  second  partial  derivatives.  In  addition  it  combines  the  desirable  features 
of  projection  methods,  conjugate  gradient  methods,  and  methods  that  solve  LP  problems 
to  obtain  feasible  directions.  Computational  results  on  a  wide  variety  of  test 
problems  are  given.  Some  comments  on  the  efficiency  of  the  algorithm  as  compared  to 
other  algorithms  is  included. 
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SUMMARY  AND  CONCLUSIONS 


Problem 


Increasingly,  the  research  community — particularly  that  inhabited  by 
management  scientists  and  operations  researchers — have  become  sensitive 
to  the  need  for  new  or  more  realistic  approaches  to  the  solution  of 
large,  complex  decision  problems.  This  new  interest  shifted  the  focus 
of  attention  from  operational,  tightly -cons trained,  single-objective 
decisions  at  lower  levels  of  organizations  to  planning  and  policy 
decisions  at  the  upper  level.  The  shift  has  not  yet  been  rewarded  with 
the  same  kind  of  dramatic  success  that  was  associated  with  highly- 
structured  operational  problems  (e.g.,  the  development  of  linear  pro¬ 
gramming).  Substantively,  the  impetus  for  this  research  was  a  problem  in 
planning  promotions  for  naval  enlisted  personnel.  Upon  investigation, 
it  became  clear  that  the  Navy’s  enlisted  advancement  system  sought  to 
achieve  multiple,  conflicting  objectives  associated  with  authorized 
strength  levels,  meeting  manpower  requirements,  minimizing  cost,  and  maxi¬ 
mizing  promotion  opportunity,  among  others.  In  attempting  to  develop  an 
optimization  technique  to  handle  such  problems,  it  was  determined  that  an 
effective  nonlinear  programming  algorithm  would  be  required.  As  a  result, 
it  was  necessary  to  address  this  problem  before  progress  was  possible  in 
the  case  of  the  enlisted  promotion  planning  problem.  Specifically, 
this  report  documents  the  development  of  a  general  algorithm  (and  the 
theory  that  underlies  it)  for  solving  the  class  of  nonlinear  programming 
problems  that  have  linear  constraints.  The  linear  constraints  can  be 
either  equations  or  inequalities  and  the  variables  can  be  free  or  non¬ 
negative.  The  objective  function  is  assumed  to  be  continuously  differ¬ 
entiable.  This  class  of  problems,  i.e.,  those  having  linear  constraints 
and  a  nonlinear  objective,  seems  to  be  particularly  appropriate  for 
manpower  management . 

Background 

Under  the  direction  of  the  Chief  of  Naval  Personnel,  this  Laboratory  is 
conducting  a  research  program  in  the  area  of  enlisted  personnel  planning. 
The  thrust  of  this  program  is  toward  the  development  of  computer-assisted 
decision  systems  for  more  effective  personnel  planning.  In  the  course 
of  this  research  it  became  apparent  that  the  ability  to  achieve  personnel 
management  objectives  is  heavily  dependent  on  actions  taken  in  the  area 
of  enlisted  promotions.  This  awareness  generated  an  investigation  of 
the  processes  underlying  the  enlisted  promotion  system,  resulting  in  the 
development  of  new  techniques  for  planning  promotions.  Still  unsolved, 
however,  was  the  problem  of  achieving  multiple  personnel  management 
objectives  which  involved  drastic  tradeoffs  in  planning  the  number  of 
promotions  in  each  pay  grade  of  each  rating  over  time.  Existing  tech¬ 
niques  available  for  the  solution  of  such  problems,  particularly  those 
of  a  non-linear  character,  indicated  some  technical  and  logical  diffi¬ 
culties  in  application. 


iii 


Approach 


During  the  last  10  to  15  years  many  algorithms  have  been  developed  for 
solving  nonlinear  programming  problems.  These  algorithms  range  from 
those  that  apply  to  special  classes  of  problems  to  some  quite  general 
algorithms  applicable  to  broad  classes  of  problems.  The  result  of  this 
development  is  a  variety  of  methods,  each  with  its  particular  strengths 
and  weaknesses.  It  is  expected  that  this  trend  in  algorithm  development 
will  continue.  However,  while  it  is  probable  that  no  universal  method 
will  emerge  that  is  clearly  superior  to  all  other  methods  for  all  prob¬ 
lems,  it  seems  quite  reasonable  to  assume  that  new  ideas  and  techniques 
will  increase  the  variety  of  methods  and  algorithms  available.  Certain 
general  techniques  have  been  used  in  solving  nonlinear  programming 
problems  and  have  been  successful  in  practice.  These  techniques 
involve  the  following  approaches:  (1)  solve  a  local  linear  programming 
problem  for  a  feasible  direction;  (2)  projection,  i.e.,  project  an 
infeasible  direction  onto  the  feasible  set;  and  (3)  solve  a  local 
conjugate  gradient  problem  for  curvature  information.  It  is  interesting 
to  note  that  algorithms  are  usually  developed  in  terms  of  only  one  of 
these  techniques  rather  than  a  combination.  However,  there  is  a  method 
that  combines  techniques  (2)  and  (3) .  The  approach  of  the  method 
developed  herein,  called  the  primal-dual  method,  is  to  combine  the 
desirable  features  of  each  of  these  techniques  into  a  single  algorithm. 
Thus,  a  large  number  of  existing  algorithms  become  special  cases  of  the 
primal -dual  algorithm.  In  addition  to  providing  a  partial  unification 
of  nonlinear  programming  algorithms,  a  very  efficient  algorithm  results. 

Findings,  Conclusions,  and  Recommendations 

The  primal-dual  algorithm  is  an  "effective"  second-order  method  in  that 
slow  convergence  is  eliminated  without  requiring  second  partial  deriva¬ 
tives.  It  has  additional  advantages  in  that  (1)  zig-zagging  due  to 
highly  eccentric  ellipsoids  is  eliminated,  (2)  if  the  stationary  point 
occurs  at  an  interior  point,  no  linear  programming  problems  are  solved 
and  the  solution  is  identified  immediately,  and  (3)  accurate  bounded 
variable  constraints  are  generated  from  the  conjugate  gradient  solution. 
Computational  results  on  a  wide  variety  of  small  (in  size)  test  problems 
indicate  the  validity  of  these  findings.  It  is  concluded  that  the 
algorithms  of  Rosen,  Goldfarb,  and  a  special  case  of  Graves  algorithm* 
can  be  viewed  as  special  cases  of  the  primal-dual  method.  A  future 
report  will  describe  the  application  of  this  method  to  the  enlisted 
promotion  problem. 


*See  References  [3]  through  [6]  and  [12]. 
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A  FIRST-ORDER,  PRIMAL-DUAL  METHOD  FOR  MINIMIZING  A 


REAL-VALUED  FUNCTION  SUBJECT  TO  LINEAR  CONSTRAINTS 


(1) 


where  is  the  number  of  equations  and  is  the  number  of  free  vari¬ 

ables.  The  function  -f ( is  assumed  to  be  of  class  C**,  i.e.,  contin¬ 
uously  differentiable.*  The  algorithm  is  a  first  order  method  in  that 
only  the  first  partial  derivatives  of  the  objective  function  are  used. 

The  conjugate  gradient  method  is  used  to  obtain  an  unconstrained 
stationary  point  when  one  exists.  The  particular  conjugate  gradient 
method  used  is  quite  efficient  in  that  no  linear  searches  are  required. 

In  the  vicinity  of  an  unconstrained  minimum  the  second-order  terms 
in  the  Taylor  series  expansion  of  *?W)  dominate.  However  in  the  vicinity 


^vectors  are  denoted  by  underlined  lower  case  letters. 


I.  INTRODUCTION 


The  purpose  of  this  paper  is  to  develop  a  general  algorithm  for 
solving  the  problem 
subject  to: 

T) 

I  %  y;  -  ri 

L  'si  ri 
yj 1  ° 


minimize : 


l  -  1 - rr\1 

L  =  i. , _ m 

j  *  ---  ^ 


f  <■  n,--- y-n) 
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of  a  constrained  minimum  the  first-order  terms  dominate.  Thus  the  method 


should  have  quadratic  convergence  (without  requiring  the  second  partial 
derivatives)  combined  with  the  ability  to  obtain  feasible  directions. 

The  algorithm  presented  here  has  the  benefits  of  a  second-order  method 
in  eliminating  slow  convergence  without  requiring  the  second  partial 
derivatives.  It  also  has  additional  advantages:  (1)  zig-zagging  due 
to  highly  eccentric  ellipsoids  is  eliminated;  (2)  if  the  stationary 
point  occurs  at  an  interior  point,  no  linear  programming  (LP)  problems 
are  solved  and  the  solution  is  identified  immediately;  and  (3)  accurate 
bounded  variable  constraints  are  generated  from  the  conjugate  gradient 
solution. 

In  order  to  present  a  complete  theoretical  development  of  the 
algorithm  we  begin  by  discussing  unconstrained  minimization  and  the 
conjugate  gradient  method.  This  is  followed  by  the  development  of  a 
first  order  accelerated  conjugate  gradient  method.  Next  it  is  shown 
that  a  feasible  direction  for  the  constrained  problem  can  be  obtained 
by  solving  an  LP  problem  and  the  method  for  determining  the  distance 
moved  is  given.  Further,  it  is  shown  that  the  Kuhn-Tucker  conditions  are 
satisfied  at  the  constrained  minimum.  A  proof  of  convergence  concludes  the 
theoretical  section.  In  the  numerical  section  a  number  of  test  problems 
are  presented.  The  algorithm  is  compared  with  existing  algorithms  by 
solving  the  test  problems  using  total  computation  time  as  the  basis  of 
comparison. 
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II.  UNCONSTRAINED  MINIMIZATION 


It  is  not  difficult  [8]  to  show  the  following. 

Theorem  1:  For  feC  let  X*  be  a  point  for  which  V-P(x*)=0*  then  "Fix) 
assumes  a  relative  minimum  at  X*  if  the  quadratic  form  (X-X*/MIX*)  (X“X*) 
is  positive  definite,or  equivalently, if  every  eigenvalue  of  HCX*J  is 
positive. 


Here  H  -  l|  j  |\  is  the  "nxn  Hessian  matrix.  It  is  also  true  that  the 

quadric  (X-.X*)  j^X-X)=-l  is  an  ellipsoid  if  all  the  eigenvalues  of 

HlxJ  r)  are  positive.  Consider  the  problem 


r\q  iivj  i  M  i  -z.  e. 


•Fix) 


x  e  E 


If  ‘ffeC^'and  X  is  an  estimate  of  X*  the  relative  minimum  then  by  Taylors 
Theorem  , 

-f  (X-)  =  f  (X*}  +  k  C**)  ^  +  r  ^  x*  , 

'A'  p 

where  ^^X -  X  .  When  the  error  term  is  zero,  minimizing  +  t  is 

equivalent  to  minimizing  the  quadratic  form  ^  HCx^^which  is  an  ellipsoid 

A” 

since  all  the  eigenvalues  of  HCx  )  are  positive.  If  X  *  X  i  is  an  estimate 

„  *  , 
of  X  then 


(xL-/)  mx+)(x-x*)  =  kL 


—  i  I  Q- , - - 
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Thus  we  want  to  find  a  sequence  of  points  such  that  the  sequence 

Geometrically  this  is  equivalent  to  finding  the  center  of  the 

ellipsoid 

We  can  extend  these  ideas  as  follows:  let  be  an  estimate  of  x* 
with-feCZ.  Then  by  Taylors  Theorem 

-fix)  = -f  V-P  (><£)•£  +  ri^^)  | 

where  ^  =  X  -  .  Let 

Fix)  = -fcxO  Vfix i)^  +  I  ^7I-I  IXt)  . 

We  now  show  that  F(x)  is  an  approximating  ellipsoid  at  Xf  If 

ftxO  -  ‘P  <-*i>  -  V-Pcxo  xL  -v  i  x'L  H  ix^  Xi  , 

then 

Fix')  =  !lxt>  4.  £  y-flx^  -  x*  M  ix^VJ  x  +  L  x '  W  (. x i')  x  . 

Now  for  any 

F  LX)  -  ?  -i  +[v-f(.x^  -X£  1-UXiW  (Xl)|  X 

+  i  (x- • 
z 

Since  a  necessary  condition  for  W(Xt)  to  be  positive  definite  is  that 

det  ^  o  »  it  follows  that  |-|  exists  [11].  Suppose  we  choose  « 

such  that 

ix^  -  x'L  McxO  +  sl'WiXi.')  =  o  . 
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Then 


®f  =  *i  ~  tVf  <  xj  } 

and 

Fcx'k  =  ?(.xL)  -  i  ot'M  ixL)  <XL  -t- 

is  an  approximating  ellipsoid  at  X^.  Since  bK^i)  is  positive  definite 

(  X  -<*/  M  (.  xiu  X  - e? )  s  O  , 

with  equality  holding  whenX-  Thus  the  minimum  of  F(>0  occurs  at 

X  *  and  is 

MiM  F(.X\  =  flX^}  ~  i  <*?  M  (.  x  L)  «  . 

2. 

The  following  has  been  established* 

Theorem  2:  If  HtX^is  positive  definite  then  the  approximating  ellipsoid 

at  *L  , 

Fix'*  =  ? C  x-L)  +  V-f  (.  xL^  -  x^  m  (.x^  j  x  +-  i  x/  H  (x^)  x 

has  the  minimum  value 

-i  , 

at  the  center 

*L+t  =  Xi.  ~  FUx  LX-J  . 
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Theorem  2  is  an  alternate  derivation  of  Newtons  method  which  finds  the 
center  of  the  approximating  ellipsoid  in  one  step.  At  this  point  we 
present  a  version  of  the  conjugate  gradient  (eg)  method  which  finds  the 
center  of  the  approximating  ellipsoid  inn  steps.  However  the  eg  method 
does  not  require  the  Hessian  or  its  inverse.  The  algorithm  given  here 
was  developed  by  Hestenes  [10].*  Its  numerical  behavior  on  the  standard 
test  functions  has  also  been  investigated  [8]. 

A  description  of  the  conjugate  gradient  method  follows. 

Starting  routine:  Select  an  initial  point  and  compute 

=  -v-Pcxj  ,  ^=3.^  . 

Main  routine:  Given  and  p  choose 

-  — K) Jk  Fk 

*  K+±  =  +  pK  ,  (  CTK  ^  O  #  arb.lrarj)  ) 

and  compute 


) 

■3k  "  ( $1“  §-k+i')/( 

.  C* 

(•K  ) 

K+l  “  —  K  +  ^  K  p  K 

.  3 

K+L  » 

> 

K 

•Pk+1 

=  +  Ph.  • 

The  basic  idea  of  accelerating  the  conjugate  gradient  method**  stems 
from  the  fact  that  the  ray  connecting  the  centers  of  two  successive 


*The  conjugate  gradient  method  was  developed  by  Hestenes  in  1952  [9]. 

**This  idea  was  suggested  to  me  by  Professor  M.  Hestenes,  UCLA. 
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approximating  ellipsoids  generally  is  a  good  direction  toward  the  minimum. 
To  develop  this  method  let  X*  be  a  relative  minimum  and  let  }  X^  and 
be  successive  estimates  of  x*.  Consider  the  Taylors  expansion  of  Xz 
about  ,  where  ^ 

-Pix^  =  -PciO  +  v-P(x1^  +  . 

>  * 

We  also  require  the  Taylors  expansion  of  X  about 

•f  tx*)  =  f  +■  V-P  i  +ii'M  i  +  r  t.  x1(  x  } 

where  £  If  ^  is  obtained  from  X^by  the  eg  method  then  by 

Theorem  2 

^  =  - Hixt'>‘v-Pc*iy  # 


We  can  now  prove  the  following  > 

,  _C 

Theorem  3:  Let  X  be  the  minimum  of  a  convex  function  tcC  .  If  XQ>  X^ } 

and  X*.  are  successive  estimates  of  X*  where  X^  is  obtained  from  X^  by  the 
eg  method  then 


■f(X+)  ■>  -fcxO  +  V-f  (  ^ 
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is  a  lower  bound  of  "f  and  if  )£.  *  £  >  O  , 

X  =  ><2.  +  j.  ^ 

II  ^ll2- 

is  the  minimum  of  -f  L  X  )  along  the  line  X  a 

Proof:  If  -f  is  convex  then  by  standard  properties  of  convex  functions 

-fix*)  2  -f  (x.j'i  +  1  ^  £tx  j  +  V? (.  ^  +  V-Pcx^  £ 

=  -PlXL)  -  V-PlfcO  M  Cx^"1  V-PiX^  +  VftX^  2 

If  Xjl.^  X  then 

tXi  + A  2)  -  -fc  X^')  <Q  >  , 

X 

and 

I  /  hA  ■PCxi  +  \2')  --PiXO  X>o  , 

X-o  ' - x - 

or 

-  V-Fcx^  2^0. 


-  i. 
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If  ^-^>Oand  ©  is  the  angle  between  and  2  then  COSQ^Oand  *f  can  be 

A  A 

reduced  along  the  line  X^X^+k^  .  Let  X  =  +  k  be  such  that 

% *  <•  * )  -  o  . 

Then 

k  =  VJ.  >  O  , 

since  ^*5^0  . 


This  theorem  is  of  no  direct  use  computationally  since  1L  is  unknown. 
However  if  the  acceleration  step  is  done  by  a  linear  search  along  ^ 
then  we  can  restrict  the  search  to  positive  values  of  k.  Limited 
computational  experience  has  shown  that  when  acceleration  can  be 
achieved  it  is  usually  quite  effective.* 


^Acceleration  works  very  well  on  Powells  function  of  4  variables 

i 

but  was  completely  ineffective  on  Rosenbrocks  parabolic  valley  [8]. 
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III.  DETERMINATION  OF  A  FEASIBLE  DIRECTION 


In  the  constrained  problem  (l)}let  X  be  the  unconstrained  stationary 
point  if  one  exists.  Let  \£°  (feasible)  be  an  estimate  of  the  constrained 
minimum,  and  A^  a  direction  emanating  from  ^°.  Thus  we  are  interested 
in  points  of  the  form  +  A  where  A^.  is  a  feasible  direction.  Now 

^°is  either  on  the  boundary  of  the  feasible  set  or  is  an  interior  point. 
To  obtain  A^_ when  is  on  the  boundary  consider  the  Taylor*s  expansion 
of  the  objective  function  } 

+  V-Pc^°)  A^i  r(  ^ .  ca) 

We  also  require  the  distance  function 

d  Ly)  =  il  x*-^.  \[z  -  ( x*-  +  —  *  C  x*  -  i 

and  its  Taylors  expansion 

d  -  <d(^°)  +  Vok^0)  +  r(  , 

where 

Vdi^°>  =  -e  t  x*-^V  • 

This  suggests  we  consider  the  convex  combination 

(.  i-  <=xg)  V-P  -  2  «B  CxA-^°)  ,  (3} 

where  O  £  ^c0  £  i  .  This  leads  to  a  stepwise  procedure  where  the  direction 
Av  is  obtained  by  solving  the  linear  programming  problem 
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subject  to: 


minimize 


L  a:;  y;  =  ^ 
I  %'  Vi  4  r‘ 

J- i 

1  V,  ‘  8* 

Vi  -  ° 


Vi 


a  -0. 


i.  =  t _ 

i=mi+l< _ tnrv 


---  ^ 


(4) 


In  order  to  assure  convergence  to  a  stationary  point  we  require  O  . 

The  upper  bound  is  found  by  choosing  so  that 

-n 

4-1 

^  (the  first  feasible  solution)  and  ^-x*.  The 
lower  bound  8^  is  found  similarly.  Note  that  for  any  solution  ^  of  the 
LP  problem  A^_is  a  feasible  direction  since  the  feasible  set  is  convex. 

The  LP  algorithm  used  here  is  the  primal-dual  algorithm  of  Graves  [5].* 
The  possible  terminal  states  of  the  LP  problems  are: 

(1)  Finite  optimal  solution  ^  -  ^-^°)  ; 

(2)  Inconsistent  constraint  (terminate  because  problem  is  inconsistent);  or 

(3)  Unbounded  solution  (does  not  occur  due  to  bounded  variable  constraints) . 


is  satisfied  when 


*It  is  particularly  efficient  in  that  no  slack  variables  or  artificial 
variables  are  required. 
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At  the  constrained  minimum  none  of  the  bounded  variable  constraints  can 
hold  as  equalities.  If  this  condition  occurs  we  simply  increase  8a.  or 
B ^  and  use  the  current  feasible  point  as  a  new  starting  point. 

In  practice  °(0is  decreased  to  zero  quite  rapidly.  When  and 

^°is  on  the  boundary  zig-zagging  can  still  arise  when  solving  (4) . 

This  is  due  to  the  eccentricity  of  the  approximating  ellipsoids.  When 
tnis  occurs  we  solve  the  local  conjugate  gradient  problem  to  estimate 
the  curvature  of  in  the  vicinity  of  ^°.  The  direction  obtained  is 

then  projected  onto  the  sub space  determined  by  the  binding  hyperplanes 
to  yield  a  feasible  direction.  The  technique  of  projection  is  similar 
to  that  developed  by  Rosen  [12].  However  instead  of  projecting  the 
gradient  we  project  one  of  the  local  conjugate  gradient  directions 
(see  the  conjugate  gradient  method) .  To  develop  these  ideas  let  B  be 
the  set  of  binding  hyperplanes  at  i.e., 

|  arui  . 

Let  B  be  a  matrix  of  rank  W  whose  rows  are  a  set  of  W  linearly  independent 
hyperplanes  from  B  .  Thus  B  is  an  ln*M  matrix  of  rank  H  .  The  subspace 
S  generated  by  the  columns  of  B  is  the  same  as  that  generated  by  the 
rows  of  0  . 

Now 

E"  =  S  +  CHS )  , 

where  CHS)  is  the  subspace  which  is  the  orthogonal  complement  of  S.  So, 
any  vector  in  E  can  be  written  uniquely  as  the  sum  of  a  vector  in  S  and 
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a  vector  in  O^S)  .  If  d  is  the  direction  obtained  from  the  local  con¬ 
jugate  gradient  problem  then 


d  -  A^P  +■  uj  j 


where  A^Pfc  OLS)  and  uj  G  S  .  Now 


giving 


^  *  Y.  xj  b-j  -  , 

S- i 


Since  A^P£  O(S)  ; 


Since  ( B  0 )  is  a  nonsingular  matrix  J 


*  *  (  b'b  f1  e'a 


giving 


l*j  -  B  ( 


and 


A^p  -  I  -  BlB'efV]  d  ^  Pd 

where 

P  =  I  -  B  (  B'B)-1^  . 
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The  matrix  P  is  called  a  projection  matrix  since  it  projects  4-  into 
OtS).  In  practice  we  calculate 

rK=M;^rK-.vft^°,PgK|r^o^  , 

.  A  A  A 

and  take  where  corresponds  to  irK  .  Actually,  any  such 

that  O  and 

vf^°)E^K^0 

A 

yields  a  direction  in  which  to  move.  It  is  not  apriori  evident  that  _pK 
is  tne  best  direction  in  which  to  move.  Now  since vftyv  , 
if  we  take  cL  -  jp  ^  and  if  Eft*  O  then 

V-P^°>Ayp  =  -d'A^p  =  -A^p/(  to)  =  -AypA^p  «.  o  . 

A 

Tnus  the  set  is  not  empty  and  V^K  exists. 

Note  that  the  calculation  of  the  projection  matrix  B  requires  a 
matrix  inversion.  In  practice  this  is  done  internally  by  the  LP 
algorithm  and  is  quite  efficient. 

If  is  an  interior  point  then  A  y  is  locally  unrestricted.  In 

this  case  no  LP  problem  is  required.  To  motivate  the  expression  for 
Ay  when  y°  is  an  interior  point  note  that  when  <=*3=0  (3)  becomes 

V- P<.^°)  A ^  . 

If  A^  -  -  then 

t^°j  Ay  -  -  U  v-Piy0)!!2”^  °  • 
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When  }  (3)  becomes 

-  2C  X*- ^°)/A^  . 

If  A^  -  (.  X+-  '^°  )  then 

=  -2iU*-^°lt^O  . 

This  suggests  we  take 

=  -  C  4-  o<x  (  x*-^°  )  ,  IS) 

where  O  ^  °<x  -  1-  .  Substituting  (5)  into  (3)  we  obtain 

-  C  i-  -o'bH  l-®<rUl  ^7-f  C^.0)llX  -  2  o<Q  o< -j-  I|  x*-  l|z 

+  i-o<8^  +  2  <*B  Vf  -  C4.) 

The  first  two  terms  of  (6)  are  negative  and  the  sign  of  the  third  term 
depends  on  the  angle  between  V-Pc^°)  and  .  Thus  in  order 

to  assure  a  gain  from  an  interior  point  when  °<q  =  0  we  require  °<x  O  . 
With  °<i  - 1  we  can  eliminate  slow  convergence  due  to  highly  eccentric 
ellipsoids.*  When  the  unconstrained  stationary  point  is  not  identified 
as  optimal  on  the  first  step  then  the  constrained  minimum  occurs  on  the 
boundary.  So  from  an  interior  point  we  want  to  move  back  onto  the  boundary 
In  practice  we  set  -  1  and  return  to  the  boundary  if  a  decrease 

in  -P<.Y)  is  achieved  on  the  boundary.  Otherwise  we  decrease  until 
°CX  -  O  .  If  can  not  be  decreased  on  the  boundary  for  °<x  -  O  then 

*A  good  example  of  this  is  given  in  the  test  problems. 
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linear  search  is  performed  along  If  a  problem  has  multiple 

local  minima  and  some  are  feasible  and  some  are  infeasible  then  an 
interesting  possibility  occurs.  Suppose  X*  occurs  at  one  of  the  infea¬ 
sible  local  minima.  If  ^is  close  enough  to  one  of  the  feasible  local 
minima  then  the  algorithm  switches  from  trying  to  find  a  point  that 
reduces  the  distance  to  X*  and  finds  the  feasible  local  minima. 

We  now  show  that  the  Kuhn-Tucker  conditions  are  satisfied  at  the 
constrained  minimum  >£*.  As  an  auxiliary  result  the  classical  Lagrange 
multipliers  are  given  in  the  final  LP  tableau.  To  show  this  observe 
that  the  Lagrangian  function  of  (1)  is 


+  X.  M  ri -  L N  Yi'l  > 

L*i  U  j  =  i  J 

where  the  A^are  called  Lagrange  multipliers.  The  necessary  conditions 
for  to  be  a  constrained  minimum  (Kuhn-Tucker  conditions)  are  (see 
Hadley  [7]  chap ter  6) 


and 


-rn 


I  x*  0.;.  - 

,  J-  i—  "i 

l  j 

*  yi 

t  ^ 

1  x*  ay  . 

j  j  =  V'i+  i>-T' 

^  t  J 

L*L  J 

Xt  *  o 

,  *■  *  vr^+i  , - rn 

(8) 

*m 

i*  i’>  - 

is**]  =°  • 

A  ’1 

19) 

(io) 


When  -  O 


and 


we  have  solved  the  LP  problem 


16 


subject  to: 


*  -  r‘ 

yj  >o 


minimize : 


Vf  <.**)£  . 


The  dual  of 
subject  to: 


this  problem  is 


yr\ 


X  * Li  XL  *  bfLy*) 
X  a\j'  X,;  i  ^  j ( y*J 


L=  L 


X-  f:  O 


maximize : 


X 


U  1, - m± 


L  -  rnL+l _ -rn 


j  =■  ^c+i.  --- ^ 


I 


-  m 


J  =  -  ^ 

t  =  i } _ rr\ 


The  necessary  conditions  (7)  and  (8)  are  satisfied  by  the  dual  con¬ 
straints  and  we  see  that  the  dual  variables  are  the  Lagrange 
multipliers.  Condition  (9)  says  that  when  an  inequality  constraint  does 
not  hold  exactly  its  corresponding  Lagrange  multiplier  *  O  .  But 
this  is  equivalent  to  the  fact  that  the  dual  variable  when  its 

corresponding  primal  constraint  does  not  hold  exactly.  Condition  (10) 
can  be  rewritten 
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Since  =  O  when  constraint  i  does  not  hold  exactly  we  have 

r*  m 

I  y?  -  L  r;  , 

j-i  A  Yj  i.«l 

or 

Vf  ly+)  y*  -  r  X*  . 

But  from  the  dual  theorem  of  linear  programming 

and  since  ^  *  )£°  +  A^  =  +  we  have 

V-F +  V-Fi^*)A^  =  r  X*  . 

It  will  be  shown  later  that  the  condition  for  termination  of  the 
algorithm  is 

I  V-f  A^  |  *  e  • 

From  these  considerations  it  is  clear  why  the  algorithm  is  called  a 
primal-dual  method.  The  feasible  direction  is  determined  from  the  primal 
problem,  and  the  identification  of  the  constrained  minimum  (by  the  Kuhn- 
Tucker  conditions)  follows  from  the  dual  problem.  In  this  framework 
Wolfe’s  algorithm  for  quadratic  programming  would  be  considered  a  dual 
method. 
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IV.  DETERMINATION  OF  THE  DISTANCE  MOVED 


We  now  determine  the  distance  moved  along  the  direction  A  ^  .  Con¬ 
sider  the  Taylor  series  expansion  of  constraint  i  about  the  feasible 
point  j 

+  Vguj°)  A  3.  , 


or 

a;¥  -  g;^°  *  SU  A*  *  ^  . 


Let 


Then  - 


A^;  =  k  A  U.  UJ  l^ere  /\  u.  -  A 

U  A^.11 


avxol.  k  y  o  . 


kaL  A sx  £r--a^°  . 

o  »  o 

Since  ^  is  feasible  ^  ^  O  .  Thus  we  have  five  cases  to 

consider  for  determining  the  restrictions  on  k  : 


1 

ri  -  =  o  }  ac  Aa  =  o 

k  unrestricted 

-  o.L^°  =  o,  S-i.  *  o 

k  unrestricted 

3 

r c  "  >  o,  i;Au. -o 

k  unrestricted 

4 

aL  Au,^  o 

k  unrestricted 

5 

L  >  o,  a^Aoi  >  o 

k  *  r;-8.l7a,A4 
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These  cases  can  be  represented  schematically  as  : 


In  case  5  let* 

A 

k  —  M  i  *4  -  S l 

1  SC  Au 

A 

Thus  the  restriction  k  £  k  ensures  us  that  the  point  ^  *•  +  ^ 

remains  feasible.  We  actually  determine  k  by  a  search  along  Au  .  The 
search  algorithm  given  here  is  also  used  for  the  acceleration  step  in  the 
eg  method. 


*The  implicit  constraints  ~  S  ;v  £  O  must  be  considered  in  calculating 
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Search  Algorithm  — 

A 

Starting  routine:  Set  (j!)  =  W  t  rv-  1  and  calculate 

+  Au  }  S  -  f  c  i1)  -  . 

If  $  *-  O  set  "n  =  2.  and  go  to  the  general  routine  . 

If  S  -  O  set  ^  "  H  l.1-  ||  . 

If  (p  —  £  set  _r'  =  i-  and  continue  . 

If  ^6  stop  k  =  O  . 

Gene ral  rout ine : 

Caluclate  +  (  k  -  Au  }and  -  f  (.  Z1”)  --(■(. 

If  $  O  set  T\  -  "n+  i.  and  continue  . 

For  S  -  O  : 

If  “  II  *  £  set  rul, 

and  continue. 

If  ^  *  II  2°-  ^.°l|  ^  f  and  ■P  t  ■g.T'  *")  -p  (. 

,  Tl  —  1 

then  V  ~ 
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V.  CONVERGENCE  OF  THE  ALGORITHM 

When  <*B-0  and  none  of  the  bounded  variable  constraints  hold 
exactly  the  LP  problem  (4)  is  equivalent  to  the  primal  problem 
subject  to:  ^ 


I  ay  Ay;  ■  rc-a-^” 

J-i 

J  =  l 


C= 

C  -  "m  j_+  1 1 _ 

j  =-nL  +  i,...-rv 


minimize : 


The  dual  of  this  problem  is 
subject  to: 

vr* 

T- 


u-l 


&  f  C^° ) 


rn 

Z  <*;:  X  •  -  -  K(*°) 

i-L  J  J  v  T 


X;.  *  O 


*J  4  O 


i  * 

i  =  t _ m 

j  =  ,  ---  r\ 


maximize: 


Z.  ^--aL^)X.  +  Z  x° 

i  = 


t*  L 


From  (2)  and  the  dual  theorem  of  linear  programming,  i.e., 


v-Fi 


vn  VI 

^°)A^  =  +  Z  y/zj 

u-1-  j  = 
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we  have 


j  -  i 


Now 


I  »!*■)*.  *  I  V°e. 


since 


.-Yn 


and 


Global  or  local  convexity  of  -ft^)  implies  that  ^1^°,  in  (11). 
Thus  the  condition  for  a  constrained  stationary  point  and  termination  of 
the  algorithm  is 


Note  that  this  condition  was  also  required  by  the  Kuhn-Tucker  conditions. 
This  leads  to  the  following 

Theorem  4:  (Convergence  Theorem)*  If  <=<0-0  then  a  stationary  point  is 
located  in  a  finite  number  of  steps. 


*This  theorem  is  a  special  case  of  a  more  general  convergence  theorem 
due  to  Graves  [6]. 


23 


Proof:  If  (12)  is  not  satisfied } then  by  the  Approximation  Theorem 

(Buck  [2]  page  243) 


•Pc^.°+k  Awt)  =  f  ly°)  +  V-Fty°)  k  All  +  r(^°  l<  Aa)  , 

where 

I  itv\  V  (  ^°f  k  Au.  )  -  O  }  u  Nil  for  hA  . 

k  4u-ro  k 

Thus  for  any  there  exists  a  $  >0  such  that  for  k  II  Aull  -  k  - 


r(*°%  k  Au.) 

From  (12)  and  (13) 


-£  k  •  €. 

2.  II  A^ll 


-  €vS 

2  II 


•f(v)  =  -P(v°)  +  Jc.  Vf(.Y0)Av  +  r  ( k  A  u.*) 


*  fi^°)  -e_^$  +  e_£  =  -  ej 

iia^ii  ziiAnju 

So  if  (12)  is  not  satisfied  then  we  can  reduce  Tlv)  by  at  least 
at  each  step  until  (12)  is  satisfied.  Thus  a  stationary  point  is 
located  in  a  finite  number  of  steps. 


113) 


•Z/zwL^W 
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VI.  DISCUSSION 


It  has  been  known  for  some  time  that  if  y_°  (an  estimate  of  the  con¬ 
strained  minimum)  is  not  close  to  y*  then  no  local  method  can  ever  be 
devised  which  obtains  a  good  direction  from  to  In  light  of  this 

fact  it  seems  advisable  to  give  up  the  search  for  better  local  directions 
and  adopt  a  radically  different  criterion,  such  as  that  embraced  in  the 
primal-dual  method.  In  the  initial  steps  we  focus  on  reducing  the  distance 
from  the  starting  point  to  the  constrained  minimum.  The  direction  from 
to  the  unconstrained  minimum  serves  as  the  surrogate  direction.  Upon 
reaching  the  boundary  we  hope  to  be  close  enough  to  the  constrained 
minimum  so  that  by  solving  (4)  with  we  obtain  .  Clearly  problems 

can  be  constructed  where  this  approach  does  not  give  a  point  on  the 
boundary  close  to  .  But  it  is  not  apriori  evident  that  we  have  moved 
further  "out  of  the  way"  than  we  might  by  some  other  method.  Although 
we  may  reach  the  boundary  at  a  point  quite  close  to  zig-zagging  can 

still  arise  due  to  the  eccentricity  of  the  approximating  ellipsoids  of 

.  Here  it  is  evident  that  second  order  information  on  the  curvature 
of  is  absolutely  necessary  in  order  to  proceed  efficiently  to  . 

The  method  presented  herein  solves  the  local  conjugate  gradient  problem 
for  the  complete  set  of  conjugate  directions  and  projects  the  "best" 
direction  onto  the  subspace  determined  by  the  binding  hyperplanes.  Since 
the  calculation  of  the  projection  matrix  requires  considerable  matrix 
operations  we  take  as  many  projection  steps  as  possible,  i.e.,  until  the 
B  matrix  must  be  changed. 


25 


The  burden  of  calculation  in  the  primal-dual  method  falls  principally 
upon  the  LP  algorithm  and  the  conjugate  gradient  algorithm.  The  LP 
algorithm  solves  (4)  for  the  feasible  direction  and  inverts  the  8  B  matrix. 
The  conjugate  gradient  algorithm  solves  the  unconstrained  problem  and 
obtains  the  complete  set  of  local  conjugate  directions  for  a  projection 
step.  Neither  the  LP  algorithm  of  Graves  [5]  nor  the  conjugate  gradient 
algorithm  of  Hestenes  [10]  are  well  known  or  widely  used.  However, 
it  would  appear  that  both  these  methods  are  the  most  efficient  of  their 
class.  Thus  it  is  hoped  that  the  primal-dual  method  presented  herein 
will  prove  to  be  the  most  effective  general  algorithm  for  non-linear 
programming  problems  with  linear  constraints.  It  is  Effectively*1  a 
second-order  method  without  requiring  second  partial  derivatives.  In 
addition  it  combines  the  desirable  features  of  projection  methods,  con¬ 
jugate  gradient  methods  and  methods  that  solve  LP  problems  to  obtain 
feasible  directions . 

In  comparing  the  primal-dual  method  with  other  methods  we  will  con¬ 
sider  (1)  a  special  case  (that  of  linear  constraints)  of  Graves  general 
nonlinear  algorithm,  (2)  Goldfarb*s  conjugate  gradient  method,  and 
(3)  Rosen*s  gradient  projection  method  for  linear  constraints.  The 
algorithm  of  Graves  is  a  first  order  method  when  only  the  first  partial 
derivatives  are  used.  The  approach  of  Graves  [6]  is  to  attempt  to  solve 
all  problems  in  this  manner.  When  slow  convergence  is  observed  in  the 
first  order  algorithm  the  second  order  conditions  are  appended  as  explicit 
constraints  provided,  of  course,  the  second  partial  derivatives  can  be 
obtained.  This  approach  expands  the  size  of  the  problem  and  requires  the 
second  partial  derivatives. 
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The  conjugate  gradient  method  of  Goldfarb  is  effective  on  problems 
with  eccentricity  in  the  approximating  ellipsoids  of  •f(^)  .  An  additional 
advantage,  as  indicated  by  Goldfarb  [3],  is  that  a  minimum  of  function 
evaluations  are  required.  This  is  important  in  some  problems  in  optimal 
control.  However,  many  problems  are  difficult  to  solve  just  due  to  size  alone 
(many  constraints  and/or  many  variables)  even  without  additional  compli¬ 
cations.  Since  Goldfarb’s  method  requires  even  more  matrix  operations 
than  that  of  Rosen  it  may  be  quite  slow  due  to  the  considerable  amount 
of  matrix  operations  on  large  matrices.  In  addition  there  is  a  large 
class  of  problems  where  simply  decreases  (unconstrained)  to-oo  ,i#e., 

there  is  no  finite  unconstrained  minimum.  When  the  conjugate  gradient 
method  is  applied  to  a  function  of  this  class  it  can  become  unstable 
quite  rapidly.  Since  Goldfarb fs  method  requires  the  conjugate  directions 
it  may  also  become  unstable  for  these  functions.  In  the  primal-dual 
method  if  instability  is  observed  when  solving  for  the  unconstrained 
minimum  we  set  and  o(x  =  o  permanently.  We  then  solve  (4)  to 

reach  the  constrained  minimum  .  If  slow  convergence  is  observed  and 
fc^pis  convex  a  direction  can  be  obtained  by  solving  the  local  conjugate 
gradient  problem.  Experience  has  shown  that  this  procedure  works  quite 
well. 

The  gradient  projection  method  of  Rosen  [12]  is  not  as  efficient  as 
Goldfarb’s  method  (see  the  numerical  results  in  [3]).  In  addition  it  is 
probably  not  as  efficient  as  Graves1  algorithm  due  to  the  amount  of  matrix 
operations  required . 
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VII.  NUMERICAL  RESULTS 

The  task  of  comparing  nonlinear  programming  algorithms  is  at 
best  a  difficult  job.  For  a  variety  of  test  problems  it  would  seem 
that  total  computation  time  provides  a  universal  norm  for  comparison. 
However,  the  complexities  of  modem  computer  operating  systems  do  not 
allow  the  same  test  problem  to  be  run  twice  with  the  same  running  time. 
Thus  computation  time  is  a  random  variable  and  the  probability  density 
function  is  unknown.  Comparing  computation  times  reported  on  different 
computers  is  even  more  difficult.  For  example,  the  test  problems 
solved  by  Goldfarb  [3]  in  1968  were  run  on  a  CDC  6600 — the  CDC  6600 
being  approximately  ten  times  faster  (depending  on  the  operating 
system  and  other  factors)  than  the  IBM  System  360  Model  65  used  here. 

The  test  problems  solved  by  Graves  [6]  in  1966  were  run  on  an  IBM  7044/ 
7094  direct  couple  system.  He  reports  that  the  time  sharing  feature  of 
that  system  could  increase  total  time  by  as  much  as  a  factor  of  four. 

His  times  include  loading  and  system  times  as  well  as  execution  time. 

Comparing  algorithms  on  the  basis  of  total  iterations  (steps)  is 
equally  unsatisfactory.  For  example  the  Goldfarb  algorithm  generally 
solves  a  problem  in  a  small  number  of  steps  although  each  step  requires 
considerable  calculation.  On  the  other  hand  the  Graves  algorithm  can 
take  far  more  steps  in  the  same  amount  of  time  since  each  step  requires 
fewer  calculations .  In  the  primal-dual  algorithm  the  situation  is  much 
more  complicated  since  a  number  of  different  kinds  of  steps  are  possible. 
In  solving  for  the  unconstrained  minimum  the  conjugate  gradient  steps 
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give  an  approximate  second  order  Newton  step,  i.e., 

Ay  =  -HCy°)'1'y-Pcy°)/  . 

In  solving  for  the  constrained  minimum  we  have 

(1)  LP  steps: 

Ay-y-y°, 

(2)  interior  steps: 

Ay  =  -  Cy0/  ,  and 

(3)  projection  steps: 

A^p  =  P  i  . 

Obviously  interior  steps  require  the  least  amount  of  calculation  and  pro¬ 
jection  steps  the  most.  The  art  of  constructing  an  efficient  algorithm 
is  in  providing  the  proper  balance  among  the  different  kinds  of  steps. 

The  algorithms  of  Graves,  Goldfarb,  and  Rosen  can  be  classified  in  terms 
of  the  type  of  step  used.  For  example,  the  Graves’  algorithm  takes  LP 
steps  with  °<6-  o  and  interior  steps  with  -  o  -  Rosen’s  algorithm 
takes  projection  steps  and  interior  steps  with  oc^  =  o  In  the  projection 
step  Rosen  projects  the  gradient,  which  is  always  the  first  conjugate 
direction  calculated.  Recall  that  in  the  primal-dual  method  the  complete 
set  of  conjugate  directions  are  calculated  and  the  "best"  one  projected. 
Goldfarb ’s  algorithm  takes  projection  steps  and  eg  steps.  In  the  projection 
step  Goldfarb  apparently  projects  the  first  conjugate  direction  which  gives 
a  gain,  i.e.,  P  pK  O  .  So  it  is  clear  that  each  of  these 

algorithms  can  be  obtained  as  a  special  case  of  the  primal-dual  algorithm. 

If  we  choose  not  to  solve  for  the  unconstrained  minimum  the  primal-dual 
algorithm  becomes  a  local  method  and  we  have  an  alternate  version  (primal- 
dual  II)  .  At  a  boundary  point  we  set  a  o  and  take  either  LP  steps  or 
projection  steps.  From  an  interior  point  we  take  conjugate  gradient  steps. 
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For  a  wide  variety  of  problems  (both  large  and  small)  the  two  primal -dual 
methods  are  probably  the  same  in  overall  effectiveness.  Of  course  for  any 
particular  problem  one  method  is  more  efficient  than  the  other.  These 
classifications  can  be  conveniently  summarized  in  the  following  table: 


Method  or 
Algorithm 

Type  of  Step* 

Boundary  Point 

Interior  Point 

Graves 

0 

9 

(D 

* 

Z)  O<x-o 

Rosen 

3)  4  =  -vft^°)/=p1 

^T-O 

Goldfarb 

3)  d.  vikere 

V-f  1^.°)  P  £K  <  O 

4) 

Primal- 

Dual 

1) 

OY~  3)  4-  ~  Csee  paje  4.40 

2.) 

Primal- 
Dual  II 

1)  ^8  =  0 

°r  3)  d-  {2|<  4  see.  pa^e  4-4) 

4) 

*1)  LP  step,  2)  interior  step,  3)  projection  step,  and  4)  eg  step. 
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The  test  problems  given  in  the  next  section  have  been  devised  to 
include  most  of  the  undesirable  phenomena  that  arise  in  constrained  mini¬ 
mization  with  linear  constraints.  In  solving  these  problems  we  will 

\ 

report  total  computation  time  and  total  iterations  (steps).  The  eight  test 
problems  are : 

(1)  constrained  ellipse  with  e-^1, 

(2)  Rosenbrock’s  parabolic  valley  with  a  linear  constraint, 

(3)  Powell’s  function  of  4  variables  with  linear  constraints, 

(4)  Fletcher  &  Powell’s  helical  valley  with  a  range  restricted  variable, 

(5)  Colville  problem  // 1, 

(6)  Colville  problem  // 2, 

(7)  Chemical  equilibrium,  and 

(8)  constrained  hyperbola  with  a  saddle  point. 

The  solution  of  test  problems  (5)  and  (6)  by  the  Graves  algorithm  as  a 
first  order  method  are  reported  in  [6].  The  solution  of  test  problems 
(5),  (6),  and  (7)  by  the  Goldfarb  algorithm  are  reported  in  [3]. 
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VIII.  TEST  PROBLEMS 


(1)  .  Constrained  ellipse  with  e~^l. 


subject  to: 


Y<- 

~Y 


Y*.  -i° 
Y* .  Y*-  -  ° 


minimize : 


ftsp  *  (.Xi-lO)2-  +  C^-5)1 


i  ooo 


start  at  (1,1). 


The  eccentricity,  e,  of  an  ellipse  is  always  between  0  and  1.  As 
e  — ^  0  the  ellipse  approaches  a  circle  and  as  e  — 1  the  ellipse  approaches 
a  straight  line.  In  this  case  the  value  of  e  is  0.9995.  This  problem 
is  designed  to  create  zig-zagging  in  ordinary  first-order  gradient  methods. 

(2).  Rosenbrock’s  parabolic  valley  with  a  linear  constraint, 
subject  to: 


Y*-  +  Y*  6  1,9 


minimize : 


start  at  (-1,2,1)  . 
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The  gradient  and  the  Hessian  are 


Vfcy)  *  (-4oo(y2-y£)  yi  -  acl-y^)  aoocy.-y*))  . 

E(3Yi-^i+°°5)  -2yi 
-EyL  1 


-  2.00 


Now  HiV  is  singular  when  column  1  is  a  linear  combination  of 
column  2,  i.e. , 

+  -005)  =  4yJ  , 


or 


yi.  ■=  -oos 


Vi*  V  +  -005 

x^a.i)  ,  fi^)x0 


^  -  tf 


Along  the  curve  y^-y^  -  .005  the  eigenvalues  of  Wty}  are  given 
by  (  Hly}  -  XI  )y.  -  Q  or 


4^-X 

~Z^ 

L-\ 

~  O 


giving 
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and  \^=  O } -  1  +  >  O  .  Since  =  0  the  approximating  ellipsoids 

2 . 

along  the  curve  ~  ^2.  “  '’•00  5  degenerate  to  parallel  straight  lines. 
This  problem  is  designed  to  foul  a  second  order  method  which  uses  the 
inverse  Hessian  and  to  create  zig-zagging  in  a  first  order  gradient 
method. 

(3).  Powell 1 s  function  of  4  variables  with  linear  constraints 
subject  to: 

yL  4  20  i  *  l,-  + 

minimize : 

tip  =  iy±  +  L0ypz  +  5t-y3-y^f  +  +  iotyi-y4}4 

start  at  (3,-l,0,l)  . 

The  gradient  and  the  Hessian  are 

=  (  2.(^1  +  loy^  +  4oCyi-y4)3,  soc^i+iONjj,  )  +  4 

Lot  ~  8(^-2  y3^  -iOty3-y4)  “4o(yt-y4)3)  , 


^il;2  +  120  Cy£-^4f, 

=  ) 

^*2  -  200  +12  Cva  -2v3)Z) 

=  O  . 

^33  =  io  +48  (y»-2y,)x| 

^4  =  *  -i2Q  lyi-y4f  , 

^44  =  io  +I20(sfl-y4f  > 

^23  =  ^32  =  ~24  , 

^14  "  ®  > 

^34  ~  ^>43  ~  * 
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The  unconstrained  minimum  of  Powell’s  function  is  obviously  zero  at 


(0,0, 0,0)  and  at  this  point 


2 


20 


0 


0 


20  200 


0 


0 


HCQ)  = 


0 


0  10  -10 


0 


0  -10 


10 


Thus  M(O)  has  rank  two  and  is  singular.  The  unconstrained  problem  is 
designed  to  create  slow  convergence  and  this  effect  carries  over  to 
the  constrained  problem. 

(4).  Fletcher  and  Powell's  helical  valley  with  a  range  restricted  variable, 
subject  to: 


O.Z  *  4  0.8 


minimize: 


loo[  (y3-ioe)z  +  ir-l'f]  y3  , 


where  y^  =  v  c.o  S  2.  TT  ©  ano(  , 

start  at  (-1,0,0)  . 

Since  y*.  r  cos  2TT  0  a  ’  ‘  “  "7  we  have  on  squaring  and  adding 
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Also  > 


UJ  Vte  r  e 


Since 


we  have 


Similarly 


and 


Since  u  = 

v*°  • 


and 


©  -  A.  +an 
2Tf 


if  s  if  dt  0  i  U.  +  <ff  i  r 

i^i.  60  clu.  i^i  ir 


if 

i^i 


if  =.  lOQO  (-X3  ~  lo  Q  )  Nfz.  4  ZOO  (  V-  1  ^  \fj 


if 


”  ^-000  C\3~10Qj  4-  200  (.r-  t  )  ^ 

ffll*l*)*]* 


if  =  200(^3-100)  •*-  2.y3  . 

^y3 


ya/yi  >  appropriate  restrictions  on  Q  are  required  so  that 

If  v  ^  o  then  Yi=0  when  G  =-4.  .  f,  3.  so  that  we  require 
J  4  >  4  4 


1  40  c  3  or  H  ^  2  IT  0  ^  OTT 

4  4’  2  "a  * 


- .k  ^  ©  <  1  or  -TT  <  e-nr  G  . 

4  4  *  2.  2- 
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n  1  "  1 


(5).  Colville  problem  // 1. 
subject  to: 

5 

ui,..  io 

J  =  i 

Yj  *  O 

minimize: 

5  5  S  S 

=  XclYj  *  X  X  Cij ;  y;y  •  *  X  ^  X 

J-l-  «-ai  jat  j=£ 

start  at  (0,0, 0,0,1)  . 

The  coefficients  are  given  as 
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1 


2 


3 


4 


5 


^Li 


be 


1 

-16. 

2. 

0.0 

1. 

0.0 

1 

o 

• 

2 

0.0 

-2. 

0.0 

0.4 

2. 

-  2. 

3 

-  3.5 

0.0 

2. 

0.0 

0.0 

-  0.25 

4 

0.0 

-2. 

0.0 

-4. 

-1. 

-  4. 

5 

0.0 

-9. 

-2. 

1. 

-2.8 

-  4. 

6 

2. 

0.0 

-4. 

0.0 

0.0 

-  1. 

7 

-1. 

-1. 

-1. 

-1. 

-1. 

1 

o 

8 

-1. 

-2. 

-3. 

-2. 

-1. 

-60. 

9 

1. 

2. 

3. 

4. 

5. 

5. 

10 

1. 

1. 

1. 

1. 

1. 

1. 

This  problem  is  designed  to  create  considerable  zig-zagging  on 
the  boundary  of  the  constraint  set  in  ordinary  first  order  gradient 
methods . 

(6).  Colville  problem  #2. 
subject  to: 

-  io  ^  ±  io  ia  4 


minimize : 

Uy)  -  Loo(yz-yl)2  *  4-  9oCy4-^3f  +  (±-y3)z 

+  +(y4-lf]  +  19-ec 

start  at  . 
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This  problem  is  an  example  of  the  case  where  the  unconstrained 
minimum  is  an  interior  point  and  hence  is  the  optimal  solution.  When 
started  from  the  given  starting  point  the  solution  path  leads  through 
(the  four-space  equivalent  of)  a  narrow,  steep-sided,  curving  valley. 

In  addition,  a  true  stationary  point  lies  in  the  path.  The  combination 
of  these  effects  is  to  create  zig-zagging  near  the  stationary  point 
which  leads  to  instability. 

(7) .  Chemical  equilibrium, 
subject  to: 


+  2.  \fa.  +  Z  y3  +  sju  ♦  Yio  =  2 


y*  +  2  y5  +  y-  +  y* 


- 1 


y3  +  yi  +  y8  +2-y=  +  y*°  -1 


yL  &  o  i*  i . io 


minimize : 


io 


Lo 


start  at  (0.1,  0.35,  0.5,  0.1,  0.35,  0.1,  0.1,  0.1,  0.1,  0.1)  . 
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The  coefficients  are  : 


C 


i 


1  -  6.089 

2  -17.164 

3  -34.054 

4  -  5.914 

5  -24.721 


6  -14.986 

7  -24.100 

8  -10.708 

9  -26.662 

10  -22.179 


The  gradient  is  easily  calculated  to  be 

4£  =.  c  •  +  I oo  /  Vi  \  i«  1, -  -  io 

**  H  ly,1 

A  detailed  discussion  of  the  physical  aspects  of  the  chemical  equilib¬ 
rium  problem  is  given  in  [1]  pp.  46-49-  Since  log  0  is  undefined  a 
numerical  problem  arises  in  the  primal-dual  algorithm.  The  constraints 
i*i, -..lO  are  in  reality  o  but  LP  algorithms 

cannot  handle  the  strict  inequality.  To  circumvent  this  problem  we  add 
the  auxiliary  constraints 

“ yi  *  - ^  l»  i, . .  to 


This  problem  is  an  example  of  the  case  where  there  is  no  finite  uncon¬ 
strained  minimum  yet  f  is  convex. 
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(8).  Constrained  hyperbola  with  a  saddle  point. 


subject  to: 

4  9 

0.5  Yl  *  V*  ^4 

yi.Vi- 2  ° 

minimize : 

=  -yty*. 


start  at  (1,1) 

This  problem  was  designed  to  trap  the  primal-dual  method  on  the  saddle 
point  at  (0,0).  It  is  not  difficult  to  show  that  the  constrained  minimum 
occurs  at  (  3/2. ,  )  and  •f  (  3/z  ,  )  -  “ZV8  .  Now  the  contours  of 

are  hyperbolas  and  the  conjugate  gradient  methods  finds  the 
center  of  hyperboloids  as  well  as  ellipsoids.  The  obvious  solution  is 
to  add  the  constraint 


-yt-'jz 


*  -e  . 
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IX.  SOLUTION  OF  THE  TEST  PROBLEMS 


(1) .  Constrained  e 

Steps 

2 

5 

1 

0 

Total  8 


llipse  with  e— ►1. 
Tz£e 

eg  steps 
LP  steps 
interior  steps 
projection  steps 


Running  Time  (sec) 

0.2 


0.5 


0.7 


Value 

unconstrained  minimum  (U.M.)  (10,  5)  0.0 

constrained  minimum  (C.M.)  (8,  5)  .004 


Lagrange  multipliers 

1  -0.4  x  10 

2  0 

3  0 


Comments:  The  conjugate  gradient  method  solved  the  unconstrained  problem 
in  1  step.  The  2nd  step  was  for  verification.  Three  of  the  five  LP 
steps  were  required  to  drive  o  •  Generally,  the  number  of  steps 

required  to  drive  is  given  by  min(m,n)  .  Each  test  problem  is 

considered  to  be  solved  when  V-Fc^°J  A 
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(2).  Rosenbrock's  parabolic  valley  with  a  linear  constraint. 


Steps 

20 

5 

1 

0 

Total  26 


Type 

eg  steps 
LP  steps 
interior  steps 
projection  steps 


Running  Time  (.sec) 

1.7 


0.5 


2.2 


U.M.  (1,  1) 

C.M.  (.96632704,  .93367296) 


Value 

0.0 

.1135190  x  10 


Lagrange  multipliers 
1  -0.2299606  x  10_1 

Comments:  The  conjugate  gradient  method  has  a  tendency  at  times  to 
take  too  large  a  step  and  become  unstable.  To  counteract  this  we  have 
introduced  an  under-relaxation  technique  [8].  Of  the  20  steps  required 
to  solve  the  unconstrained  problem,  4  included  under-relaxation. 
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(3).  Powell's  function  of  4  variables  with  linear  constraints. 


Steps 

6 

3 

6 

0 

31 

Total  46 


Type  Running  Time  (sec) 

eg  steps 

acceleration  steps  1.1 

LP  steps 

interior  steps  10.3 

projection  steps  _ 

11.4 


Value 

U.M.  (0,  0,  0,  0)  0.0 

C.M.  (.50332384,  -.045560064,  .23582560,  .30641062)  0.11378385 


Lagrange  multipliers 


1  -  .4010324 

2  0. 

3  0. 

4  0. 

5  0. 


Comments:  Since  the  projection  steps  were  consecutive  only  one  projection 

matrix  was  calculated,  thereby  making  these  steps  quite  efficient. 
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(4).  Fletcher  &  Powell's  helical  valley  with  a  range  restricted  variable. 


Steps 

Type  Running  Time  (sec) 

16 

eg  steps 

3 

acceleration  steps  2.4 

5 

LP  steps 

1 

interior  steps  0.3 

0 

projection  steps 

Total  25 

2.7 

U.M.  (1,  0,  0) 

Value 

0.0 

C.M.  (0.8,  0,  0) 

4.0 

Lagrange  multipliers 

1  -40. 

2  0. 
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(5)  Colville  problem  #1. 


Steps 

6 

8 

0 

2 

Total  16 


Type  Running  Time  (sec) 

eg  steps  1.2 

LP  steps 

interior  steps  2.6 

projection  steps  _ 

3.8 


Value 

U.M.  (.5242668,  .8826745,  1.258478,  .7411507,  .3355687)  -61.44833 

C.M.  (0.3,  0.3334676,  0.4,  .42831015,  .22396490)  -32.348679 


Lagrange  multipliers 


1 

0. 

6 

-11.839545 

2 

0. 

7 

0. 

3 

-5.1740399 

8 

0. 

4 

0. 

9 

-  .10389623 

5 

-3.0611088 

10 

0. 

Comments:  This  problem  was  solved  by  both  Graves  and  Goldfarb.  Graves* 
algorithm,  as  a  first  order  method,  did  not  do  well  because  of  zig-zagging. 
It  would  appear  that  the  primal-dual  method  solved  this  problem  roughly 
twice  as  fast  as  Goldfarb* s  algorithm. 
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(6).  Colville  problem  #2. 


Steps 

Type 

Running  Time  (sec) 

Total  41 

eg  steps 

6.6 

Value 

U.M.  (C.M.)  (1,  1,  1,  1)  0.0 

Comments:  This  problem  was  also  solved  by  both  Graves  and  Goldf arb . 
Graves’  algorithm,  as  a  first  order  method,  was  trapped  at  the  non- 
optimal  stationary  point.  Goldfarb's  algorithm  reached  the  global 
minimum  in  105  steps.  Although  it  was  considerably  slow  in  the  vicinity 
of  the  stationary  point  it  remained  stable.  The  primal-dual  algorithm 
became  unstable  at  the  14th  step,  in  the  vicinity  of  the  stationary 
point,  However,  this  step  resulted  in  movement  to  a  point  from  which 
the  global  minimum  was  reached  in  27  steps. 

The  results  of  this  test  problem  should  not  be  taken  too  seriously. 
The  conjugate  gradient  method  has  the  desirable  property  of  moving  away 
from  true  stationary  points,  but  if  forced  to  calculate  in  the  immediate 
vicinity,  it  tends  to  become  unstable. 
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(7) .  Chemical  equilibrium, 


Steps 

^  1 
11 
0 
54 

Total  65 


Type  Running  Time  (sec) 

eg  step  0.3 

LP  steps 

interior  steps  57.4 

projection  steps  _ 

57.7 


No  finite  unconstrained  minimum 


Best  Known 
C.M.  Solution  [13] 


Value 


1  .040587 

2  .146323 

3  .785045 

4  .001409 

5  .485317 

6  .000692 

7  .027265 

8  .017867 

9  .036919 

10  .095984 

-47.761079 


1  .040668 

2  .147730 

3  .783153 

4  .001414 

5  .485247 

6  .000693 

7  .027399 

8  .017947 

9  .037314 

10  .096872 

-47.761106 
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Lagrange  multipliers 


1  -  9.789526 

2  -12 .97198 

3  -15.22707 


Comments:  This  is  the  only  test  problem  of  the  eight  that  provided  any 
particular  difficulty  for  the  primal-dual  method.  The  difficulty 
probably  arises  from  the  fact  that  although  is  convex,  it  decreases 

(unconstrained)  to  -  ©o  .  The  best  known  solution  could  have  been 
achieved  if  a  smaller  zero  level  had  been  chosen.  It  is  interesting  that 
Goldfarbfs  algorithm  solved  this  problem  with  relative  ease  (about  8  times 
faster  than  the  primal-dual  algorithm) .  It  is  likely  that  the  Goldfarb 
algorithm  generates  different  conjugate  directions  than  that  given  here.* 


*However,  theoretical  work  for  a  quadratic  function,  indicates  that, 
Hestenes  and  Fletcher-Powell  both  generate  the  same  set  of 
conjugate  directions. 
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(8).  Constrained  Hyperbola  with  a  Saddle  Point. 


Steps 

1 

4 

0 

0 

Total  5 


Type  Running  Time  (sec) 

eg  step  0.1 

LP  steps 

interior  steps  0.1 

projection  steps  _ 

0.2 


saddle  point  (-.470924  x 
constrained  saddle  point 


_11  _11 
10  ,  -.470924  x  10  ) 

(0,  0) 


Value 

0. 

0. 
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In  terms  of  some  of  the  standards  described  by  Ignizio,*  some  additional 
information  on  the  numerical  part  of  this  study  is  warranted;  namely: 

Date  of  computation — July  1972  to  October  1972 
Computer  utilized — IBM  System  360 /Model  65 
Programming  language — FORTRAN  IV 

Total  number  of  problems  attempted — In  addition  to  the  test  problems 
reported,  a  number  of  other  small  problems  were  solved.  Since  they 
exhibited  no  interesting  mathematical  properties  and  were  solved  quite 
rapidly  they  were  not  reported.  Each  test  problem  is  considered  to  be 
solved  when  V-f(.^°)  . 

Computation  time — The  running  times  were  calculated  by  the  FORTRAN 
external  function  ITIME.  The  timer  was  set  up  to  be  the  first  and  last 
executable  statements  in  the  program.  The  times  reported  are  CPU  times 
only . 


*Ignizio,  J.,  "On  the  Establishment  of  Standards  for  Comparing  Algorithm 
Performance,"  TIMS  Interfaces,  2(1),  Nov.  1971,  pp.  8-11. 


51 


REFERENCES 


[1]  Bracken,  J.  and  McCormick,  G. ,  Selected  Applications  of  Nonlinear 
Programming  (John  Wiley  &  Sons,  New  York,  1968). 

[2]  Buck,  R.,  Advanced  Calculus  (McGraw  Hill,  New  York,  2nd  Ed.,  1965). 

[3]  Goldfarb,  D.  and  Lapidus,  L.,  "Conjugate  Gradient  Method  for 
Nonlinear  Programming  Problems  with  Linear  Constraints," 

I&EC  Fundamentals,  Vol.  7,  142-151  (1968). 

[4]  Goldfarb,  D. ,  "Extension  of  Davidon's  Variable  Metric  Method  to 
Maximization  Under  Linear  Inequality  and  Equality  Constraints," 

SIAM  J.  Appl .  Math.,  Vol.  17,  No.  4  (July  1969). 

[5]  Graves,  G.,  "A  Complete  Constructive  Algorithm  for  the  General 
Mixed  Linear  Programming  Problem,"  Naval  Res.  Logistics  Quart., 

Vol.  12,  No.  1,  1-34  (March  1965). 

[6]  Graves,  G.  and  Whinston,  A.,  "The  Application  of  a  Nonlinear 
Programming  Algorithm  to  a  Second  Order  Representation  of  the 
Problem,"  Cahier  de  Centre  d1 Etudes  de  Recherche  Operationelle 
Bruxelles,  Belgique,  Vol.  11,  No.  2  (1969). 

[7]  Hadley,  G. ,  Nonlinear  and  Dynamic  Programming  (Addis on -Wesley, 
Reading,  Massachusetts,  1964) . 

[ 8 J  Hatfield,  G. ,  "Unconstrained  Minimization  of  a  Real-Valued  Function," 
unpublished  report,  UCLA,  Nov.  1971. 

[9]  Hestenes,  M.  and  Stiefel,  E.,  "Methods  of  Conjugate  Gradients  for 
Solving  Linear  Systems,"  Journal  of  Research  of  the  National 
Bureau  of  Standards,  Vol.  49,  No.  6  (1952). 

[10]  Hestenes,  M.  ,  "Multiplier  and  Gradient  Methods,"  Journal  of 
Optimization  Theory  and  Applications,  Vol.  4,  No.  5,  303-320 
(Nov.  1969). 


53 


[11]  Noble,  B.,  Applied  Linear  Algebra  (Prentice-Hall,  1969). 

[12]  Rosen,  J.,  "The  Gradient  Projection  Method  for  Nonlinear  Program¬ 
ming,  Part  1.  Linear  Constraints,"  SIAM  J.  Appl.  Math.,  Vol.  8, 
181-217  (1960). 

[13]  White,  W. ,  Johnson,  S.,  and  Dantzig,  G. ,  "Chemical  Equilibrium  in 
Complex  Mixtures,"  J.  Chem.  Phys.,  Vol.  28,  751-755  (1958). 


54 


If 


DISTRIBUTION  LIST 


Chief  of  Naval  Operations  (OP-01) 

(OP-09 8T) 

(OP-099) 

(OP-39) 

(OP-964) 

Commandant  of  the  Marine  Corps,  G-l  Division,  Code  A01M 
Chief  of  Naval  Training  (Code  N-33) 

Chief  of  Naval  Material  (MAT  03P) 

Chief  of  Naval  Personnel  (Pers-llb) 

(Pers-A12)  (5) 

(Pers-A3)  (25) 

Commander  Antisubmarine  Warfare  Force,  U.  S.  Atlantic  Fleet  (Code  70) 
Office  of  Naval  Research  (Code  458)  (2) 

Naval  Postgraduate  School 
Naval  Examining  Center 
Center  for  Naval  Analyses 

Naval  Personnel  Research  and  Development  Laboratory 
Office  of  the  Assistant  Secretary  of  the  Navy 
(Manpower  and  Reserve  Affairs) 

Chief  of  Research  and  Development,  U.  S.  Army 

Army  Research  Institute  for  Behavioral  and  Social  Sciences 

Personnel  Research  Division,  Lackland  AFB 

Personnel  Management  Development  Office,  0P0,  U.  S •  Army 

Human  Factors  Operations  Research  Laboratory,  Bolling  AFB 

Science  and  Technology  Division,  Library  of  Congress 

Defense  Documentation  Center  (Attn:  DDC-TC)  (12) 


55 


U 1 520'  7  7 


